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Abstract 

We study the spin-spin correlation function in or near the T = ground 
state of the antiferromagnetic Ising model on a triangular lattice. At zero 
temperature its modulation on the sublattices gives rise to two Bragg peaks 
in the structure factor, and a known expression for the algebraic decay of 
correlations enables us to examine the form of the diffusive scattering. We 
do so by means of a comparison between exact results and data calculated 
using standard Monte Carlo techniques. At non-zero temperatures the finite 
correlation length alters this form, and we account for the change by proposing 
a generalisation of the zero temperature pair correlation function. The size 
dependence of our simulation data is investigated through a novel finite-size 
scaling analysis where t = e _2//T is used as the temperature parameter. 
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Since Onsager's evaluation of the free energy for the Ising-model on a rectangular lattice 
]1|] the corresponding model on a triangular lattice has served as an afterthought requiring 
some modification of the analysis |J. The Triangular Antiferromagnetic Ising (TAI) model 
has also attracted attention in its own right, it being a simple manifestation of a fully 
frustrated system ||. Very recently, the interest in the TAI model has been revived through 
its near-perfect experimental realisation in the yavapaiite layered structure of anhydrous 
alums such as RbFe(S04)2 @. 

Mathematically the isotropic TAI model to be studied in this paper is defined by the 
Hamiltonian 

H = J^SiSj, (1) 

where J > is the antiferromagnetic exchange coupling and Sj = ±1 are Ising spins on the 
triangular lattice depicted in Fig. |I[ The label indicates a sum over nearest neighbour 
pairs each pair being counted only once. Frustration in the ground state arises from the 
inability of the system to simultaneously satisfy the "local packing rule" that three spins 
on an elementary triangle are pairwise each other's nearest neighbour, and the "global 
packing constraint" that, in order to minimise energy, the ground state must have as many 
antiferromagnetically satisfied nearest neighbour bonds as possible. In other words, it is 
impossible to orient three spins in a pairwise antiparallel fashion. It can thus be shown 
that the ground state is macroscopically degenerate with a finite entropy per spin s = 
\ k B Jo /3 ln(2 cosu) dcu ~ 0.323 k B . 

In Wannier's original approach || the method of Onsager (diagonalisation of the transfer 
matrix using representation theory on an associated Lie algebra) was modified to include 
the diagonal interactions of the TAI model. Due to the technical complexity of this method 
much work was done to achieve a simplification of the algebra involved (see references in 
H), the result being the now well-known technique of reducing the problem to chains of 
interacting fermions |J. 

At the same time an alternative scheme, known as the combinatorial method, began to 
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emerge. Using topological theorems Kasteleyn ||[7j evaluated the configurational generating 
function of the problem covering a lattice with dimers. The success of this approach relies 
on the fact that the combinatorics can be lumped into a mathematical entity known as 
the Pfaffian, which is simply the square root of a determinant. The Ising problem is then 
addressed by counting dimer configurations on a modified lattice, in which every vertex has 
been substituted by a suitably oriented polygon in order to ensure a correct counting of the 
states summed over in the Ising partition function [0,|J. 

This scheme (as opposed to the algebraic method) is easily augmented to include the 
evaluation of correlation functions || expressed as Toeplitz determinants. Using asymptotic 
properties for these objects, Stephenson M provided the long-distance behaviour of the two- 



spin correlations along the three main directions for the TAI model with |IJ or without flT 
anisotropy. 

A recent interest fl2Hl8l has been taken in generalising the Hamiltonian (|l|) to include 
ferromagnetic next-nearest neighbour interactions as well as anisotropy. Such additional 
couplings tend to lift the ground state degeneracy, since, dividing the triangular lattice 
into the usual three sublattices, next-nearest neighbour pairs belong to the same of these 
sublattices (see Fig. |T]). By applying mappings which are constructed to automatically 
satisfy the antiferromagnetic nearest neighbour constraint, the ground state ensemble can 
be investigated within the context of solid-on-solid (SOS) |T^JT3|JT5|] or domain-wall |T6l . |P7 
models. 

In this paper we shall concern ourselves with the nature of the correlations in the low- 
temperature TAI model. Despite the fact that the abundance of entropy prevents the emer- 
gence of a long-range order, even at zero temperature [[|, the existence of some kind of 
short-range order or pattern formation should be evident from a visual inspection of Fig. |2|, 
and the T = state is indeed a critical one as witnessed by the algebraic decay of correlations 
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We shall study this ordering by means of the structure factor 
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S(p) = J2e ip - r (s s r ), (2) 



which is proportional to the cross section for quasi-elastic single scattering |[L9|| . When 



scattering, for example, neutrons from a rare-gas monolayer adsorbed on graphite |2(J the 
ordering is directly displayed through the positions and shapes of peaks in the structure 
factor. 

By Fourier transforming the known T = expression for the pair-correlation function 
|TT| we find good agreement with numerical data for the structure factor as calculated from 
standard Monte Carlo (MC) data. Also, simulations performed at non-zero temperature 
render unchanged values for the Bragg points but a broadening of the line shape, which 
makes us suggest a generalisation of the expression for the pair-correlation function to the 
T > case. Finally, the role of finite-size effects in the MC data is discussed by means of a 
finite-size scaling analysis where t = e~ 2//T is used as the temperature parameter. 



Stephenson [] 1 made an asymptotic expansion of the spin-spin correlation function for 



the T = TAI model 



cos f — ) 

(s s r ) ^M-forr>l, (3) 



valid along the three main directions of the lattice. Since the derivation of this result 
depends heavily on certain asymptotic properties of the Toeplitz determinants involved, a 
similarly stringent result valid for arbitrary lattice directions could hardly be worked out 
along these lines. However the symmetry between the tree sublattices suggests that the 
general expression is formed by replacing the factor cos(^p) with a weight factor having the 
value +1 when so and s r are on the same sub lattice, and — | otherwise fll3 |. 



Recently it has been demonstrated that a staggered field on the sublattices can be 
mapped onto a period-6 spin wave operator within the context of an equivalent SOS-model 
fl"3| . By renormalisation the latter is reduced to a Gaussian model thus allowing for a deter- 

(s) 1 

mination of the associated critical exponent as Xq = j which in in nice agreement with the 
algebraic decay ~ r~ 2X e in Eq. (H). Although this of course corroborates the original as- 



sumption of a staggered field we are not aware of anyone having formerly verified it through 
a direct study of the TAI model. 



where we have introduced the radially symmetric function /(r) = l/\/r for the purpose of 
later generalisation. It is seen that the role of the weight factor is to produce two Bragg 
peaks at wave vectors qi and q 2 , as shown in Fig. |3|. Interestingly, these positions are also 
predicted by a simple mean field theory |2(J , pleading no knowledge of such exact results as 
Eq. (Q). Furthermore, the knowledge of an algebraic decay makes it possible to determine 
the shape of the diffusive scattering. 

Since the expression for (sos r ) is only valid for r ^> 1 we expect our evaluation to work 
for |p — qj\ <C 7r (the extent of the first Brillouin zone). With a reasonably fast p-space 
decay of the peaks this restriction will have negligible consequences. The conversion of the 
sum to an integral will cause no trouble for a sufficiently large system, but when making 
computer simulations on modestly sized lattices we must be on guard for finite-size effects. 

We have performed standard Monte Carlo simulations on L x L lattices for different sizes 
up to L = 900. Periodic boundary conditions are imposed and we demand L to be divisible 
by 6 in order to avoid the introduction of screw dislocations in the ground state ||15|| . After 
initialising the lattice paramagnetically we equilibrate it at a temperature T (measured in 
units of ks/ J) by performing a suitable number of MC steps per spin (MCSS). 

Monitoring the excess energy per spin at T = as a function of MCSS (see Fig. |j) 
we find that the system after a short transient time enters a regime of linear relaxation 
where it looses most of its excitational energy. Qualitatively this regime corresponds to 
a 'trivial' relaxation during which the energy can be efficiently lowered by updating the 
spin configuration locally. After a certain cross-over time which, as demonstated on the 
inset of Fig. f|, is proportional to I a new regime of roughly algebraic decay is entered 



Defining qi j2 = (±-f-, 0) it is easily checked that | (e 



tqi-r _j_ e iq 2 -r^ j g wanted weight 



factor, and the structure factor is evaluated as 




(4) 
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as the annealing defects generated during the quench are slowly eliminated and long-range 
correlations are built up. Although the system at this point is energetically very close to 
the ground state, we must take the length of our runs to be considerably longer than this 
cross-over time in order to obtain reliable data for the Bragg peaks in the structure factor. 
To ensure that the MC data for different lattice sizes are of a comparable quality we shall 
thus take the number of MCSS to be proportional to L and fix MCSS at 10,000 for L = 900. 

After thermalisation we calculate the structure factor S(p) for the wave vectors in the 
vicinity of p = qi which are compatible with translational invariance, reinitialise, and make 
a new quench. (Since only a small fraction of the simulation time is spent in the regime 
of linear relaxation, as mentioned above, it would not be advantageous to perform the 
simulations as sequential heating runs. Instead, the method of repeated thermal quenches 
employed here guarantees the independence of the different runs.) For each value of T and 
L we average S(p) over 20 quenches. Interference with the finite system size, however, 
manifests itself as a kind of noise on the scale of the least allowable wave vector, even after 
calculating the ensemble average. This ruggedness tends to obscure the shape of the peak, 
wherefore we dispose of it by doing a mild coarse-graining replacing each value of S(jp) with 
the average of itself and its four nearest neighbours. The result is a smooth peak as shown 
in Fig. |. 

The position of the Bragg peaks as calculated from their first moment compares 
favourably with the theoretical T = result. For T = 0.5 and L = 900 we find 
qi = 7r (1.33340(6), —0.00010(9)) ~ (if>o) with similar results for other temperatures and 
lattice sizes. 

To avoid a circumstantial notation we focus our attention on one of the two peaks in 
S(p), picking out the j — 1 term of Eq. (|]) in the following. Shifting the p-space origin to 
the centre of the peak we have 

//*oo 
d 2 re ip '7(r) = 2ir j dr rf{r) J (pr), (5) 

where Jo(x) is a spherical Bessel function of order zero. 
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Consider first the case of T = where Stephenson's result is represented by f(r) = 
l/\/r. The integral can be written in terms of gamma functions |23| as 5'(p) ~ 
2nV2p-^T (§) /r (J) or 

5(p) ~ 3.0033 (6) 

with T] = ~. 

This prediction is readily verified from our simulation data by computing a circular 
average of S(p). The result for the L = 900 lattice at T = is displayed in Fig. g. As 
expected we find a power law behaviour S(p) ~ p-^-wo)^ valid in all but the immediate 
vicinity of the peak. The exponent is 77900 = 0.45. 

Finite-size effects enter in two different ways. Firstly, the exponent deviates slightly from 
the theoretical value rj = \. Secondly, the impossibility of reproducing exact divergences on 
a finite lattice means that the power law fit breaks down near the peak centre p = 0. This 
can be interpreted as the effect of a large but finite correlation length as described below. 

At finite temperatures T > the system is no longer critical and is, as opposed to 
the zero-temperature case, characterised by a finite correlation length. The Stephenson 
expression for the decay of correlations in this case JTT], Eq. (1.31)] can be rewritten as 

p-r/| 

f(r) = ^, (7) 



where the correlation length is given by £ 1 = — In tanh ■=, thus implying that v — 1 [21 
Actually the oscillatory factor is now proportional to cos (c(T)^p), but since c(T) — > 1 as 



T — > the deviation from cos is of no consequence for sufficiently low temperatures 
22fl . Note that Eq. (0) reduces to the T = result for £(T) — ► 00 as it should. 
Performing the r-integral [E3] we arrive at 



S(p)~ 3.0033p-( 2 -^(pO, (8) 
where the function g(x), shown in Fig. [7], can be expressed as a hypergeometric function 



24 1 . Since g{x) — > 1 for x — > 00 the consistency with Eq. (y") is evident. 



For high values of £ the expression for S(p) retains its pure power-law form for all wave 
vectors p, since those allowed are bounded from below by On the other hand we have 
g(x) ~ x 2 ~ v for x — > 0, whence a lower value of £ makes S(p) tend towards a constant for 
small wave vectors. 

Our data for the circular averages of the diffusive scattering at non-zero temperatures are 
shown in Fig. ||, and the agreement with the predictions of Eq. (|8|) is seen to be quite good. 
The curves shown are parametrised by the values of £ (T) obtained from the exact expression 
given above, but as shown in Table 1 these values are within 10 % of the optimum parameters 
found from a non-linear least-square fit. However, at very low temperatures (T < 0.25) when 
the deviation from a pure power law decay becomes minute we are unable to reproduce the 
T — > divergence in £ (T) due to the finite size of our lattices. 

We now turn our attention to a scaling analysis of the data obtained for various lattices. 



According to the standard theory of finite-size scaling p5| the spin-spin correlation function 
should scale with system size L as (s s r ) = L~ 2/3 ^ u f tL 1 ^^, where t = -^-^ is the reduced 
temperature. When applying this scaling ansatz to the TAI model two caveats must be taken 
into account. Firstly, t is not a suitable measure of the dimensionless temperature in a model 
with T c = 0. Secondly, the critical exponent (3, which is defined in terms of the spontaneous 
magnetisation below the critical temperature, is a very doubtful parameter indeed, since its 
domain of definition simply does not exist. 

We shall dispose of the first complication by appealing to the Coulomb gas description 
given in Ref. [|13|] according to which a temperature induced excitation from the TAI model 
ground state is equivalent to the formation of a vortex in the associated SOS-model. These 
vortices correspond to all three spins on an elementary triangle being aligned and accordingly 
have the Boltzmann weight e~ 2//T when expressed by our dimensionless variables. The second 



problem is easily eliminated, since an inspection of Kadanoff's block spin argument |26 
reveals that the factor L~ 2/3 ^ u is only conventional and may be replaced by L~ r] without ever 
referring to (3. 

Finally, we define a sublattice independent correlation function T(r, T) = (sos r )(35i r ^A — 
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2), where i r is the label of the sublattice (we fix the origo so that io = A), i.e., by eliminating 
the dependence on the weight factor. The finite-size scaling hypothesis now assumes the form 

T(r,T)=L-H( K j J e-^ T L 1 ^y (9) 

Note that 77 is the correlation function exponent pertinent to the staggered field. Albeit the 



scenario of Ref. |13[] opens the possibility of three diffent spin waves, each with its own value 
of 77, the period-2 wave would only be present if we were to include an external magnetic 
field, and the period-3 wave would only be relevant if the coupling constants were sublattice 
dependent. 

At zero temperature we have, from Eq. ([3D, that T(r, 0) ~ r~ v for r ^> 1, whence 

f{x,0)~ x - r ' (10) 

for values of x = £ satisfying 0<i< 1. As witnessed by Fig. |9] this asymptotic behaviour 
of f(x, 0) is in fact nicely brought out for all values of x shown on the graph, at least for the 
small lattices. When increasing the lattice size we first see an unexpected fall-off for large x 
and eventually, for L = 900, even a clear tendency for T to ramify in three individual curves 
corresponding to the three sublattices. 

The latter observation could be taken to imply that the largest lattices have been insuf- 
ficiently thermalised to correctly bring out the T = correlations for distances above some 
one hundred lattice constants. To check this hypothesis we have produced some extra runs 
for the L = 300 lattice, using a different number of MCSS than in our main series of data. 
In Fig. [H^ we have redisplayed the L = 300 data as above along with the results obtained 
when using 10%, 30% and 300% of the nominal number of MCSS. The results clearly cor- 
roborate our suspicion, that both the anomalous fall-off and the sublattice upsplitting are 
indeed effects of an insufficient thermalisation time. Considerations on currently available 
CPU time, however, prevented us from performing longer runs than the ones reported here. 

The dependence of the scaling function f(x, y) on its second variable can be examined 
by plotting correlation data for a range of non-zero temperatures and different lattice sizes 
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in a series of curves with a fixed value x = xq [27j. In Fig. |TT] we show the plot of f(xo,y) 
versus y = e~ 2//T L 1// ^, and it is found that the choice u — 1.0 ± 0.1 makes the data collapse 
on distinct curves parametrised by the different values of xq- As expected from the above 
discussion on insufficient thermalisation the quality of the Xq — | curve is inferior to that of 
the other branches of the graph. 

Note, that for sufficiently low y, i.e., at low temperatures, f(x ,y) is constant for a 
fixed value of Thus T(r, T) ~ L~ v ~ r~' n , whence the zero temperature algebraic decay 
T(r, 0) ~ r~ r] is in fact valid even for small finite values of e~ 2//T L x l v '. 

In summary, we have shown that Stephenson's expression (Eq. |3|) can indeed be applied 



to all lattice directions in the way suggested by Ref. |13j. Moreover, we have confirmed the 
validity of the appropriate generalisation to small non-zero temperatures by exhibiting the 
agreement with our Monte Carlo data. Finally, a finite-size scaling analysis demonstrated 
that the correct way of approaching the T = critical point is through the temperature 
parameter t = e~ 2 ^ T . 

The authors wish to thank Ole Mouritsen for having originally brought the TAI model 
to our attention. 
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FIGURES 

FIG. 1. The L x L triangular lattice spanned by ai,2 = (±|,^) can be divided into three 
interpenetrating sublattices labelled A, B and C. Sublattice A is spanned by af 2 = (±§, 
When imposing periodic boundary conditions, L must be divisible by six in order to preserve this 
sublattice structure and avoid the introduction of screw dislocations in the ground state. Note, that 
next-nearest neighbours belong to the same sublattice, whence a ferromagnetic coupling between 
them tends to align the spins on each sublattice and thus lift the ground state degeneracy. 

FIG. 2. MC simulation of the TAI model on a 90 x 90 lattice clearly shows the difference 
between the appearance of the paramagnetic initial state (left) and the labyrinthine patterns of 
the locally ordered ground state (right), which is reached after 1000 MC steps per spin. Here, spin 
up (down) is represented by the presence (absence) of a dot. 

FIG. 3. A point in the reciprocal lattice spanned by bi^ = 27r(±l, -^), with its six nearest 
neighbours. The inner hexagon is the first Brillouin zone. Variational mean field theory predicts 
that the Bragg points describing the continuous transition to an ordered phase be located according 
to the Lifshitz criterion, i.e., that they be points of high symmetry in the first Brillouin zone. 
Apart from the origo, which characterises a genuine long-range order, the possibilities satisfying 
this symmetry demand are marked by a dot in the figure. Imposing the additional condition that 
the Landau free energy be minimal, we are left with the two Bragg peaks at qi and q2- Fourier 
transformation of the pair-correlation function confirms the mean field scenario. 

FIG. 4. The main graph shows the excess energy per spin at T = 0, measured relative to the 
ground state, as a function of Monte Carlo time for a range of different lattice sizes L = 300, 600, 
900, 1200 and 1500. After a short transient time the system enters a regime of linear relaxation ter- 
minating in a cross-over to an algebraic decay as the ground state is approached. As demonstrated 
in the inset this cross-over time is proportional to L. 
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FIG. 5. Close-up of the structure factor in the vicinity of the Bragg point qi = ^-^,0J for a 
system of size L = 900 at temperature T = 0.5. An average over 20 independent quenches as well 
as a mild coarse-graining have been performed. 

FIG. 6. Circular average showing the line shape of the T = Bragg peak simulated on a 
L = 900 lattice. The exponent of the algebraic decay is —(2 — r/goo) = —1.55. 

FIG. 7. The function g(x) |24| controlling the deviation of the diffusive scattering from a pure 
power law form. From the asymptotic behaviour we infer the constancy of S(p) at small values of 
x = pC (using g{x) ~ x 2_r? for x — > 0) and the unchanged form at large ones (using g(x) — > 1 for 
x — » oo). 

FIG. 8. Fits to the shape of the diffusive scattering at temperatures T = 0.25, 0.35, 0.40, 
0.45, 0.50 and 0.60 as indicated by the labels. For clarity the graphs have been shifted along the 
5(p)-axis. 

FIG. 9. Finite-size scaling at T = 0. The graphs for f(x, 0) versus x for different system sizes 
(L = 90, 150, 300, 600 and 900 as labelled) collapse on a universal curve. The breakdown of scaling 
at large x is due to an insufficient thermalisation time. 

FIG. 10. Replotting the L = 300 data for different thermalisation times (the label indicates the 
percentage of the nominal number of MCSS) we infer that the deviation from the universal scaling 
function f(x,0) ~ x _r? could be ameliorated by making longer runs. 

FIG. 11. Finite-size scaling at T > 0. The graphs for f(xo,y) versus y fall on distinct branches 
according to the value of the parameter xq = r/L. For clarity the different branches have been 
shifted along the ordinate. 
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TABLES 



T 


0.25 


0.35 


0.40 


0.45 


0.50 


0.60 


e (fit) 


172 


124 


82 


48 


29 


15 


f (exact) 


1490 


152 


74 


43 


27 
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TABLE I. Correlation lengths £ (in units of the lattice constant) for a L = 900 lattice as a 
function of temperature T. 



15 



.Jlllllilir 



jIlllSlllllllllll^ 



1 




2 




-2.5 -2.0 -1.5 -1.0 -0.5 

log(p) 



\ 

IP 



100.00 



10.00 - 



0.10 - 



0.01 



-I — 



~ I — I — I I I I 





i — 


50 




T 


r — 




x 


Lw 


/ r = 


20 




- L/ 


/ r = 


15 


X 


- L/ 


/ r = 


12 




: L/ 


/ r = 


10 


o 


"L 7 


/ r = 


6 


A 



xxx 
□ 

o 



o 



o 



A 



A 



X 



A 



+ 



A a A 
A 



A 



^A 



+ 

<X + 



o 



o 



□ 

o 
o 



A 



A 



_i i 



_i i i i i_i 



0.01 



0.10 



1.00 
e- 2/T V /v 



10.00 



x 
+ 

x 

□ 
o 



100.00 



